/*
 * ConjugateGradient.h
 *
 *  Created on: Mar 24, 2011
 *      Author: curt
 */

#ifndef CONJUGATEGRADIENT_H_
#define CONJUGATEGRADIENT_H_

#include "CUDA_vec.h"
#include <cmath>
#include "utility.h"


class ConjugateGradient {
public:
	template<typename Matrix> void solve(const Matrix & A, CUDA_vec<float> & x, const CUDA_vec<float> & b, float tol)
	{
		typedef CUDA_vec<float> Vector;
		int N = x.getSize();
		float dk, dkp1, alphak, bdelta;
		Vector pk(N);
		Vector rk(N);
		Vector sk(N);
		//Initial search direction = b - A*x
		compute_residual(A, x, b, pk);
		compute_residual(A, x, b, rk);
		dk = dotProd(pk,pk);
		bdelta = dotProd(b,b);
		while(dk > tol * tol *bdelta)
		{
			A(pk, sk);
			alphak = dk / dotProd(pk, sk);
			x.AXPBY(1,alphak,pk);
			rk.AXPBY(1,-alphak,sk);
			dkp1 = dotProd(rk,rk);
			pk.AXPBY(dkp1/dk, 1, rk);
			dk = dkp1;
		}
	}
};



#endif /* CONJUGATEGRADIENT_H_ */
